iconEuler Examples

Stability of BDF

by R. Grothmann

We investigate the stability of the BDF (backward differential formula) following the work of Gear e.a.

A BDF is a formula to solve differential equation. It uses multiple steps, and it is implicit. We assume the solution in the points x_0,...,x_n given as y_0,...,y_n and determine y_{n+1} implicitely such that

Stability of BDF

where p interpolations (x_i,y_i) for i=0,...n, and satisfies

Stability of BDF

First, let us determine the formula for p for the case n=1.

>function p(x) &= a+b*x+c*x^2;

It suffices to solve the special case

Stability of BDF

since we are only interested in equidistant discretization here.

>sol &= solve([p(-h)=y0,p(0)=y1,diffat(p(x),x=h)=yd2],[a,b,c])
                     h yd2 + 2 y1 - 2 y0      h yd2 - y1 + y0
       [[a = y1, b = -------------------, c = ---------------]]
                             3 h                      2
                                                   3 h

We define a function for the solution.

>function ps(x) &= p(x) with sol[1]
         2
        x  (h yd2 - y1 + y0)   x (- h yd2 - 2 y1 + 2 y0)
        -------------------- - ------------------------- + y1
                   2                      3 h
                3 h

The value of this function at h defines p(x_{n+1}) from the previous values and the derivative p'(x_{n+1}).

>function ystep(y0,y1,yd2,h) &= ps(h)|ratsimp
                         2 h yd2 + 4 y1 - y0
                         -------------------
                                  3

Here is an Euler function, which solves differential equations using this two step BDF. We use iteration for the solution of the non-linear equation

Stability of BDF

This is not optimal and must be replaced with faster solvers in real life.

>function bdf (f,a,b,n,y0) ...
   x=linspace(a,b,n);
   y=zeros(size(x));
   y[1:2]=runge(f,x[1:2],y0);
   h=(b-a)/n;
   for i=2 to n;
   y[i+1]=y[i]+f(x[i],y[i])*h;
   repeat
     yold=y[i+1];
     yd2=f(x[i+1],yold);
     y[i+1]=&:ystep(y[i-1],y[i],yd2,h);
     until yold~=y[i+1];
   end;
   end
   return y;
 endfunction

Let us show that the order of this procedure is 2.

>error1000=bdf("-2*x*y",0,1,1000,1)[-1]-1/E
4.9057824697e-007
>error2000=bdf("-2*x*y",0,1,2000,1)[-1]-1/E
1.22635597177e-007
>log(error1000/error2000)/log(2)
2.00010545601

Next, we solve the special equation

Stability of BDF

The non-linear problem for p can be solved explicitely for this equation. We get an explicit equation for y1 given y_0,y_1 to y_2.

>difgl &= solve(ystep(y0,y1,lambda*y2,h)=y2,y2)
                                y0 - 4 y1
                        [y2 = --------------]
                              2 h lambda - 3

Let us show for lambda=-5, that this recursion does indeed converge to the correct solution.

>lambda=-5; n=1000; h=1/n; ...
 sequence("(x[n-2]-4*x[n-1])/(2*h*lambda-3)",[1,exp(lambda/n)],n+1)[-1]
0.00673766561896
>exp(-5)
0.00673794699909

For stability, a formula must be able to solve this case, even if h*lambda is large. We want the iteration to converge to 0 for each lambda<0, or even better, for each complex lambda with re(lambda)<0. This is called A-stability.

The two step PDF is stable for this lambda and h, as we see with the following start values.

>sequence("(x[n-2]-4*x[n-1])/(2*h*lambda-3)",[1,0],n+1)[-1]
-0.00340277570995
>sequence("(x[n-2]-4*x[n-1])/(2*h*lambda-3)",[0,1],n+1)[-1]
0.0101912705026

We can study the problem exactly. The sequence of y_n is defined by a so-called difference formula. To solve such formulas, we can make the ansatz y_n=t^n. It follows easily, that solutions of this formula are determined by the zeros of

Stability of BDF

If t is a zero, t^n will satisfy the difference formula. We are interested to see |t|<1, if re(lambda*h)<0.

Replacing t=1/z, we seek the solutions of

Stability of BDF

We want to show that no |z| <= 1 with re(p(z))<0 exists.

>function q(z) &= 3-4*z+z^2
                              2
                             z  - 4 z + 3

Let us plot the image of |z|=1. We see that indeed, the method is A-stable, since the interior of the image is in the right half plane.

>z=exp(I*linspace(0,2pi,1000)); plot2d(q(z),r=10):

Stability of BDF

Discussing that involves a discussion of the function

>&realpart(q(cos(x)+I*sin(x)))|trigsimp
                       2         2
                  - sin (x) + cos (x) - 4 cos(x) + 3

We can do that analytically, but here is a simple plot.

>plot2d("2*x^2-4*x+2",-1,1):

Stability of BDF

Three term BDF

Next we do the same for the three term BDF.

>function p(x) &= a+b*x+c*x^2+d*x^3;
>sol &= solve([p(-2*h)=y0,p(-h)=y1,p(0)=y2,diffat(p(x),x=h)=yd3],[a,b,c,d]);
>function ps(x) &= p(x) with sol[1];
>function ystep(y0,y1,y2,yd3,h) &= ps(h)|ratsimp
                    6 h yd3 + 18 y2 - 9 y1 + 2 y0
                    -----------------------------
                                 11

>sol &= solve(ystep(y0,y1,y2,lambda*y3,h)=y3,y3)
                           - 18 y2 + 9 y1 - 2 y0
                     [y3 = ---------------------]
                              6 h lambda - 11

The formula for

Stability of BDF

has three terms now. We follow the same ideas as above.

>function q(z) &= 11+(num(rhs(sol[1])) with [y0=z^3,y1=z^2,y2=z])
                           3      2
                      - 2 z  + 9 z  - 18 z + 11

This polynomial now stretches into the left half plane. This makes the method unstable for h*lambda=i.

>plot2d(q(z),r=50):

Stability of BDF

An instable Formula

We derive an instable and unusable formula. The idea is to interpolate

Stability of BDF

Stability of BDF

using Hermite interpolation, then to set

Stability of BDF

>function p(x) &= a+b*x+c*x^2+d*x^3
                           3      2
                        d x  + c x  + b x + a

>sol &= solve( ...
   [p(-h)=y0,diffat(p(x),x=-h)=yd0,p(0)=y1,diffat(p(x),x=0)=yd1],[a,b,c,d])
                               h (2 yd1 + yd0) - 3 y1 + 3 y0
        [[a = y1, b = yd1, c = -----------------------------, 
                                             2
                                            h
                                         h (yd1 + yd0) - 2 y1 + 2 y0
                                     d = ---------------------------]]
                                                      3
                                                     h

We can easily derive a formula for the next y-value.

>function ystep(y0,yd0,y1,yd1,h) &= p(h) with sol[1]
        h (2 yd1 + yd0) + h (yd1 + yd0) + h yd1 - 4 y1 + 5 y0

To investigate the stability, we use y'=lambda*y as above.

>y2 &= ystep(y0,lambda*y0,y1,lambda*y1,h)
        h (2 y1 lambda + y0 lambda) + h (y1 lambda + y0 lambda)
                                           + h y1 lambda - 4 y1 + 5 y0

Set h*lambda=c.

>y2c &= ratsimp(y2 with lambda=c/h)
                     (4 c - 4) y1 + (2 c + 5) y0

This difference recursion should be stable for re(c)<0. So let us investigate this.

>&c with solve(y2n-y2c,c), function h(t) &= % with [y2n=t^2,y1=t,y0=1]
                          y2n + 4 y1 - 5 y0
                          -----------------
                             4 y1 + 2 y0


                              2
                             t  + 4 t - 5
                             ------------
                               4 t + 2

We now do not want to have values |x| > 1, which touch the negative axis.

>z=exp(I*linspace(0,2pi,1000)); ...
 plot2d(h(z),r=6):

Stability of BDF

But we see that such x exist for all c with re(c)<0. Indeed, the algorithm fails completely, even for moderate c.

>lambda=-1; n=100; h=1/n; c=lambda*h
-0.01
>y=sequence("(2c+5)*x[n-1]+(4c-4)*x[n-2]",[1,exp(-lambda*h)],n+1)[-1]
-2.03142847388e+057

Examples